Libraries

library(ggplot2)
library(vegan)
library(dplyr)
library(purrr)
library(readr)
library(stringr)
library(plotly)

Utility constants and functions defined

## A set of colors
color.list <- c("#E64B35FF", "#4DBBD5FF", "#00A087FF", "#3C5488FF", "#F39B7FFF", "#8491B4FF",
                "#91D1C2FF", "#B09C85FF", "#FAFD7CFF", "#82491EFF", "#B7E4F9FF", "#FB6467FF",
                "#526E2DFF", "#E762D7FF", "#FAE48BFF", "#A6EEE6FF", "#95CC5EFF")
## ggplot2 preset theme
figtheme <- theme_classic() +
  theme(text = element_text(size=23,face='bold'),
        axis.title.y=element_text(margin=margin(0,15,0,0)),axis.title.x=element_text(margin=margin(15,0,0,0)),
        plot.margin = unit(c(1,1,1,1), "cm"),
        plot.title = element_text(margin=margin(0,0,15,0), hjust=0.5))
theme_set(figtheme)
## Auxilary function to merge two profiles
merge2 <- function(x, y){
  aux <- function(x) read_tsv(x, col_names=c("tax", str_remove(basename(x), '[\\._\\-].*')))
  if(is.character(x)){
    dat.x <- aux(x)
  }else{
    dat.x <- x
  }
  dat.y <- aux(y)
  full_join(dat.x, dat.y)  
}

Read and merge profiles

profile.list <- read_lines(params$profile_list)
tax.profile <- reduce(profile.list, merge2)
tax.profile <- tibble::column_to_rownames(tax.profile, 'tax')

tax.profile[is.na(tax.profile)|tax.profile<params$presence_absence_thre] <- 0 ## impute NAs and set low abundance taxon to 0
tax.profile <- tax.profile[rowSums(tax.profile) != 0, ] ## remove empty taxa
tax.profile

\(\beta\)-diversity overview with Bray-Curtis dissimilarity

distance <- vegdist(t(tax.profile))
cmds <- cmdscale(distance, eig = TRUE)
perc <- (cmds$eig/sum(cmds$eig))[1:2]*100
plot.dat <- data.frame(cmds$points)

p <- ggplot(plot.dat, aes(x=X1, y=X2)) + 
    geom_density_2d(color='grey') + 
    geom_point(size=3) + 
    scale_color_manual(values=color.list[c(1,4)]) + 
    labs(x=sprintf('PCoA 1 [%.1f%%]', perc[1]), y=sprintf('PCoA 2 [%.1f%%]', perc[2])) 
ggplotly(p)

Analysis session information

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: CentOS release 6.5 (Final)
## 
## Matrix products: default
## BLAS/LAPACK: /mnt/software/unstowable/miniconda3-4.6.14/envs/shotgunMetagenomics_r_v3.6.0/lib/libopenblasp-r0.3.7.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] plotly_4.9.1    stringr_1.4.0   readr_1.3.1     purrr_0.3.2    
## [5] dplyr_0.8.0.1   vegan_2.5-6     lattice_0.20-38 permute_0.9-5  
## [9] ggplot2_3.1.1  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_0.2.5  xfun_0.11         splines_3.6.1     colorspace_1.4-1 
##  [5] htmltools_0.3.6   viridisLite_0.3.0 yaml_2.2.0        mgcv_1.8-28      
##  [9] rlang_0.3.4       pillar_1.3.1      later_1.0.0       glue_1.3.1       
## [13] withr_2.1.2       plyr_1.8.4        munsell_0.5.0     gtable_0.3.0     
## [17] htmlwidgets_1.5.1 evaluate_0.13     labeling_0.3      knitr_1.26       
## [21] httpuv_1.5.2      crosstalk_1.0.0   parallel_3.6.1    Rcpp_1.0.1       
## [25] xtable_1.8-4      scales_1.0.0      promises_1.1.0    jsonlite_1.6     
## [29] mime_0.6          hms_0.4.2         digest_0.6.18     stringi_1.4.3    
## [33] shiny_1.3.2       grid_3.6.1        tools_3.6.1       magrittr_1.5     
## [37] lazyeval_0.2.2    tibble_2.1.1      cluster_2.0.8     crayon_1.3.4     
## [41] tidyr_0.8.3       pkgconfig_2.0.2   MASS_7.3-51.3     Matrix_1.2-17    
## [45] data.table_1.12.6 assertthat_0.2.1  rmarkdown_1.12    httr_1.4.0       
## [49] R6_2.4.0          nlme_3.1-139      compiler_3.6.1